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1.  Introduction 


As  discussed  in  Scientific  Report  #1  (Seitter  and  Colby  1989;  hereafter  SC), 
the  goal  of  this  research  project  is  to  create  a  meso-0  scale  prediction  model  which 
is  capab.e  of  being  run  operationally  on  a  relatively  small  computer.  The  approach 
being  taken  is  centered  on  recasting  the  equations  in  a  new  form  which  allows  the 
critical  <  ipects  of  the  atmospheric  boundary  layer  to  be  included  without  a  large 
number  o  vertical  layers  in  the  model.  This  allows  significantly  fewer  grid  points 
in  the  three-dimensional  domain  which  leads  to  increased  computational  speed.  This 
aspect,  along  with  other  consistent  approximations  in  the  equations  and 
parameterizations  provides  a  means  of  producing  a  model  which  is  many  times  faster 
than  other  mesoscale  models  currently  in  use  today  [see  Appendix  B  of  Pielke  (1984) 
for  information  on  the  computational  aspects  of  many  current  mesoscale  numerical 
models],  A  high  resolution  model  can  only  be  implemented  on  a  regional  domain 
given  the  current  limitations  of  computers.  The  influence  of  changes  in  the  larger 
scale  weather  patterns  are  included  by  nesting  the  high  resolution  model  in  a  model 
with  coarser  horizontal  resolution  (but  the  same  vertical  structure).  This  coarse 
grid  model,  in  turn,  can  receive  boundary  information  from  the  operational  models  of 
the  National  Weather  Service. 

The  basic  model  formulation,  as  originally  conceived  and  discussed 
extensively  in  SC,  consisted  of  a  two-layer  “boundary  layer  coordinate"  model  which 
was  formally  nested  in  the  lowest  layer  of  a  four-layer  <7 -coordinate  model.  While 
the  initial  testing  reported  in  SC  showed  promising  results,  additional  and  more 
critical  tests  indicated  the  vertical  nesting  procedure  could  lead  to  accumulating 
errors  that  eventually  destroyed  the  solution.  Further  work  showed  that  the 
nesting  procedure  failed  to  conserve  energy  in  the  finite  difference  equations,  and 
this  lack  of  energy  conservation  appeared  to  contribute  to  some  aspects  of  poor 
.nodel  performance.  The  model  structure  described  in  SC  was  therefore  abandoned 
in  favor  of  that  described  in  this  report.  While  direct  comparisons  with  the 
previous  model  formulation  will  not  be  discussed  here,  it  should  be  noted  that  tne 
current  model  formulation  has  outperformed  the  previous  formulation  in  every 
aspect  of  testing  and  is  only  slightly  more  costly  in  terms  of  computation  time. 


The  current  model  formulation  will  be  discussed  in  section  2,  along  with 
changes  and  additions  to  the  boundary-layer  parameterizations  discussed  in  depth  in 
SC.  Section  3  will  present  the  results  of  model  simulations  testing  the  new 
formulation  under  a  variety  of  circumstances,  and  some  sensitivity  studies  involving 
the  specification  of  the  horizontal  diffusion.  Some  aspects  of  the  prototype 
operational  formulation  for  the  model  are  discussed  in  section  4.  Comparisons  of 
test  simulations  of  the  boundary  layer  parameterization,  including  radiation 
calculations,  with  other  boundary  layer  work  will  be  presented  in  section  5,  and  the 
report  concludes  with  section  6. 


2.  Model  Description 

a.  Basic  model  equations 

Most  mesoscale  models  rely  on  a  large  number  of  layers  near  the  surface  in 
order  to  resolve  explicitly  the  growth  and  decay  of  the  planetary  boundary  layer 
and  to  correctly  simulate  the  fluxes  of  heat  and  moi'-ture  that  couple  the 
atmosphere  to  the  surface.  This  leads  to  models  with  a  large  number  of  layers  in 
the  vertical,  resulting  in  large  numbers  of  computations  per  timestep  and  relatively 
long  computation  times  even  in  supercomputer  environments  (Pielke  1984;  Appendix 
B).  On  the  other  hand,  some  early  models  showed  success  in  simulating 
nonboundary-layer-driven  flows  (such  as  mountain  lee-waves)  with  relatively  few 
layers  (Anthes  and  Warner  1978),  and  other  models  looking  specifically  at  boundary 
layer  processes  were  able  to  successfully  capture  the  important  features  by 
treating  the  boundary  layer  as  a  single  layer  which  could  dynamically  grow  and 
collapse  (Lavoie  1972;  Colby  1983). 

The  present  formulation  seeks  to  marry  these  two  approaches  into  a  single 
three-dimensional  mesoscale  model.  The  lowest  layer  of  the  model  is  the  boundary 
layer  —  it  is  allowed  to  grow  in  depth  or  collapse  at  each  grid  point  as  the 
simulation  proceeds  in  response  to  the  surface  fluxes  produced  in  the  boundary- 
layer  and  radiation  parameterizations.  The  layers  above  the  boundary  layer  adjust 
dynamically  to  the  changing  boundary-layer  depth,  while  simulating  the  horizontal 
and  vertical  advections  of  momentum  and  thermodynamic  variables  and  maintaining 
the  proper  balances  that  hold  at  the  meso-/3  scale. 
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The  vertical  coordinate  which  allows  for  this  changing  boundary-layer  depth 
is  a  modification  of  the  common  o-coordinate  [cr  =(p  —  pt)/(ps  —  pt)  where  p«  is  the 
surface  pressure  and  pt  is  a  pressure  level  specified  as  the  top  of  the  model  (we 
take  p  £  =  1 00  mb  for  this  study)].  If  we  let  ah  represent  the  height  of  the  top  of 
the  boundary  layer  (see  Fig.l),  we  can  define  a  new  vertical  coordinate,  TJ ,  as 


V  = 


cr  — 
H 


H 


o*  a  < 

1  — c rh  a  >  ah 


(2.1) 


We  refer  to  T]  as  “boundary  layer  coordinates."  The  vector  momentum  equation, 
hydrostatic  relation,  continuity  equation,  thermodynamic  equation,  and  specific 


y=-l.O 


P  =  Pt 


(H‘l?)  =  0 


Fig.  1.  Schematic  structure  of  the  five-layer  ^-coordinate  model. 
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humidity  conservation  equation  can  be  written  in  the  T]  system  as 


arffV  ,  3 utHV  ,  3utHV  ■  drrWHf] 
at  dx  dy  dT] 


—  ir HV<j>  —  irHaVo-T  —  / k  X  irHV  +  TtHFuv 


where  it  —  Ps—pu  Q  is  the  specific  humidity,  and  the  other  terms  have  their 
normal  meteorological  definitions  (see  the  Appendix  for  information  on  the 
derivation  of  these  equations  and  a  list  of  terms).  It  should  be  noted  that  the  c r 
which  appears  on  the  right-hand  side  of  (2.2)  is  a  dependent  variable  which  will  not 
be  constant  on  a  constant  T]  surface,  in  general.  Expansion  of  this  gradient  of  a 
product  results  in  an  additional  pressure-gradient-force  type  term  not  present  in 
traditional  cr-coordinate  models.  It  should  also  be  noted  that  the  humidity  is 
treated  here  as  a  passive  scalar,  so  (2.6)  simply  represents  conservation  of  q  (that 
is,  dq/dt  —  0,  except  for  the  eddy  diffusion  term  added  to  the  rhs)  in  the  rj 
system.  When  condensation  and  evaporation  are  included,  appropriate  source  and 
sink  terms  will  need  to  be  added  to  (2.6)  and  another  conservation  equation 
governing  the  condensate  will  be  required. 

Equation  (2.4)  is  not  used  directly  in  the  model,  but  it  does  provide  the 
means  of  calculating  vertical  velocities  and  the  rate  of  change  of  surface  pressure. 
Integration  of  (2.4)  from  the  model  top  to  the  surface  yields 

-  -  j>-  +  4w,’ru)]dT7 


9 7T 
3 1 


(2.7) 
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where  Ht  and  H-  represent  the  values  of  H  above  and  below  the  boundary  layer 
top,  respectively,  as  given  by  (2.1).  Integration  of  (2.4)  from  the  top  of  the  model 
down  to  a  specific  interface  level,  along  with  the  use  of  37r/3t  found  with  (2.7) 
allows  the  determination  of  (Hr?)  at  each  interface  as 

(ot,  _  -m+n  +  Iff]  -  (2.® 

The  Appendix  verifies  that  the  continuous  equations  transformed  into  77- 
coordinates  satisfy  the  classical  energy  conservation  constraints  of  the  atmosphere. 
It  is  desirable  to  have  the  finite  difference  forms  of  the  equations  satisfy  these 
energy  conservation  constraints  as  well,  and  the  Appendix  outlines  how  the  vertical 
differencing  must  be  carried  out  to  preserve  these  properties.  It  also  shows  how 
the  energy  constraints  lead  to  specific  forms  for  the  solution  of  the  thermodynamic 
equation,  (2.5),  and  the  integration  of  the  hydrostatic  equation,  (2.3),  to  find  the 
geopotentials  of  midlevels  of  the  77-layers.  The  resulting  finite  difference  forms 
are  considerably  different  in  structure  from  equations  (2.2)-(2.6),  and  also  different 
from  the  forms  introduced  in  SC. 

The  last  terms  on  the  rhs  of  (2.2),  (2.5),  and  (2.6)  include  a  “friction”  term,  F, 
which,  above  the  boundary  layer,  is  given  by  a  horizontal  eddy  diffusion.  This 
term  is  modeled  by  a  simple  Fickian  diffusion,  KV2^,  where  K  is  a  constant  eddy 
viscosity  and  <t>  is  the  variable  of  interest  (u,  v,  T,  or  <?).  In  most  simulations,  we 
let  K  —  2  x  iO6  m2  s-!  for  the  momentum  components  and  K  =  0  for  temperature 
and  moisture  (see  section  3  for  a  more  complete  discussion  of  diffusion  in  the  model 
and  the  results  of  some  sensitivity  tests). 

h.  Grid  domain  and  horizontal  nesting 

A  staggered  grid  is  used  in  both  the  vertical  and  horizontal  directions.  In 
the  vertical,  all  variables  are  layer  quantities  except  vertical  velocities,  which  are 
defined  at  interface  levels  (see  Fig.  1).  Horizontally,  velocities  are  defined  on 
staggered  points  which  surround  the  points  on  which  all  other  variables  are  defined. 
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In  order  to  increase  the  overall  model  domain  size  and  move  the  lateral  boundaries 
away  from  the  area  of  primary  interest,  a  horizontal  nesting  of  the  model  is 
employed  as  developed  by  Zhang  et  al.  (1986).  A  fine  grid  mesh  (FGM)  with  20  km 
resolution  is  nested  in  a  coarse  grid  mesh  (CGM)  with  60  km  resolution.  A  3:1  ratio 
of  FGM  points  to  CGM  points  is  necessary  with  a  staggered  grid  so  that  both 
“velocity”  and  “thermodynamic"  points  can  be  coincident  in  the  overlap  region 
(Zhang  et  al.  1986).  The  CGM  domain  covers  1320  km  x  1320  km  while  the  FGM 
domain  is  480  km  x  480  km  for  “thermodynamic”  points  (all  displays  will  be  made  on 
“thermodynamic”  point  arrays,  with  any  displayed  velocities  being  averaged  to  these 
points).  As  described  in  somewhat  more  detail  in  SC,  the  two-way  interactive 
nesting  procedure  of  Zhang  et  al.  (1986)  is  used  with  a  few  minor  modifications.  In 
the  calculation  of  tendencies  for  the  “thermodynamic”  points  in  the  FGM,  for 
example,  a  simple  linear  interpolation  is  used  between  CGM  points  nearest  the 
boundary  FGM  point  rather  than  the  “Lagrangian  interpolation”  used  by  Zhang  et 
al.  (1986). 

Zhang  et  al.  (1986)  discuss  a  Newtonian  damping  scheme  which  is  applied  near 
the  FGM/CGM  interface  to  help  control  noise  resulting  from  the  overspecification 
of  the  pressure  tendencies  there.  As  will  be  shown  in  later  sections,  we  discovered 
similar  noise  generation,  especially  on  outflow  boundaries  where  the 
overspecification  is  most  severe.  We  therefore  plan  to  include  an  additional 
smoothing  mechanism  near  the  boundaries,  but  not  in  the  form  described  by  Zhang 
et  al.  (1386).  Zhang  et  ai.  (1986)  describe  their  Newtonian  damping  scheme  as  being 
given  by  (for  a  o-coordinate  model) 

=  •••  -OrK  -  W)/Td(a)  (2.9) 

where  the  ir-weighted  velocities  on  the  rhs  are  meant  to  be  taken  at  the  “latest 
time  level,”  and  the  overbar  represents  a  two-point  average  using  values  on  either 
side  of  the  point  of  interest  (which  is  located  one  grid  interval  in  from  the  FGM 
boundary).  The  quantity  Td(cr)  is  a  relaxation  constant  which  is  specified  as  a 
function  of  height  [varying  between  20At  and  lOOAt  for  Zhang  et  al.  (1986)1.  We 
will  look  at  this  noise  control  technique  in  some  more  detail. 
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Suppose  we  consider  applying  this  Newtonian  damper  to  a  u-momentum 
variable  U  (where  U  might  be  irHu  in  ^-coordinates  or  ir u  in  a -coordinates)  which  is 
governed  by  an  equation  of  the  form 


du  =  F 


(2.10) 


Then,  without  damping,  a  value  of  U  at  timestep  n-(-l  at  gridpoint  (i,j)  would  be 
given  by  (using  the  leapfrog  scheme) 


Utj  =  Uun~l  +  2A  tFun 


(2.11) 


If  the  damping  3cheme  is  applied  to  this  point  as  implied  by  Zhang  et  al.  (1986)  and 
Kurihara  and  Bender  (1980),  we  obtain  (assuming  an  east  or  west  boundary) 


U  tJn+l  =  Uijn~l  +  2A  tFun 


2Alfrr  "+t  1 

-  tJ  -  ^ 


lu,+1Jn+‘  + 


U,-u"+')] 


(2.12) 


This  equation  is  semi-implicit,  but  can  be  applied  explicitly  as  long  as  it  is  not  also 
applied  at  points  (i  +  l,j)  or  ((—  1,J),  because  new  values  of  U i+i/+1  and  Ui_lJn+l  can 
be  found  first,  allowing  (2.12)  to  be  solved  for  UtJn+1. 

If,  instead,  we  use  (2.11)  to  produce  an  “undamped”  estimate  of  t/1Jn+1,  denoted 
L/,/"  ',  we  can  rewrite  (2.12)  as 


0,J n+l  -  ^[0un+'  -  \{Ui+un+x  +  i/i-. 


(2.13) 


It  is  easy  to  see  that  (2.13)  represents  a  simple  numerical  filter  applied  at  time  level 
n-f-l  (Shapiro  1970),  with  the  quantity  (2Af /Td)  serving  as  the  nondimensional  filter 
factor.  Since  the  same  effect  can  be  obtained  through  the  use  a  Fickian  diffusion 
term  (Kurihara  and  Bender  1 980),  we  have  chosen  to  abandon  the  Newtonian  damping 
term  in  favor  of  an  increased  diffusion  coefficient,  K,  on  those  grid  points  near  the 
boundary. 
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c.  Time  integration  and  outer  lateral  boundary  condition 

Time  integration  for  the  model  is  performed  using  the  leapfrog  scheme  with 
an  Asselin  filter.  The  time  step  for  FGM  points  is  20  s  and  for  COM  points  is  60  s. 
The  flow  relaxation  condition  of  Davies  (1976)  is  used  on  lateral  boundaries 
following  the  work  of  Seitter  (1987)  who  found  that  this  condition  was  well-behaved, 
provided  a  simple  means  of  allowing  external  information  to  be  introduced  into  the 
model,  and  did  not  require  the  smoothie-  operator  necessary  in  the  Perkey  and 
Kreitzberg  (1976)  sponge.  The  flow  relaxation  condition  requires  a  5  gridpoint  wide 
region  near  the  boundary  for  application,  and  solutions  in  this  “relaxation  region” 
should  be  considered  modified.  The  flow  relaxation  condition  is  only  applied  at  the 
lateral  boundaries  of  the  CGM  domain,  however,  so  FGM  points  are  not  significantly 
influenced  by  the  relaxation  region. 


d.  Coupling  of  the  model  with  the  boundary-layer  parametrization 

The  basic  equations  which  make  up  the  boundary-layer  parameterization  are 
described  in  detail  in  SC  and  will  not  be  reproduced  here.  The  boundary-layer 
package  provides  the  model  with  the  fluxes  of  heat,  moisture,  and  momentum  from 
the  surface,  and  calculates  the  rate  of  change  of  boundary-layer  height  which  is 
fundamental  to  the  ^-coordinate  formulation.  A  summary  of  these  quantities  is 
shown  in  Table  1. 

It  is  important  to  note  that  the  boundary-layer  routine  diagnoses  the  current 
boundary  layer  depth  (hydrostatically)  from  the  current  ah  when  it  is  called,  and 
returns  a  time  rate  of  change  of  boundary-layer  height,  dh/dt,  where  h  is  in 
geometric  height  above  the  ground.  The  model  requires  both  the  rate  of  change  of 
boundary-layer  height  and  the  new  boundary-layer  in  terms  of  ct,  so  a  conversion 
must  be  made.  The  hydrostatic  relation  may  be  written  in  cr-coordinates  as 


9  <t> 
da 


—  trot 


(2.14) 


which  can  be  integrated  from  the  surface  (a=l)  to  the  top  of  the  boundary  layer 
(cr  =<7„)  to  obtain 
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Table  I.  Quantities  provided  to  the  model  by  the  boundary-layer 
parameterization. 


Variable  returned 


Description 


dh 

dt 

dTG 

dt 

dGW 

dt 

Tkb 

Qkb 

SH 

LH 

GS 


Rate  of  change  in  height  of  boundary  layer  top 
Rate  of  change  of  ground  temperature 
Rate  of  change  of  ground  wetness 

Diagnosed  “surface"  temperature 
D’agnosed  “surface”  mixing  ratio 
Sensible  heat  flux 
Latent  heat  flux 
Surface  soil  heat  flux 

Net  radiation  (incoming  shortwave  minus 
outgoing  IR) 

Surface  stress  (friction) 


—it 


a  da 


(2.15) 


With  the  excellent  approximation  that  0  =<  gh,  and  making  use  of  the  mean  value 
theorem,  we  can  obtain 


h  =  1  -  c t„)  (2.16) 

Taking  the  partial  derivative  with  respect  to  time  of  this  expression  and 
rearranging  yields 

= _ Q_dh  ,  (1  -  credit  ,  (1  -  <rh)dd  n 

dt  Tradt  'r  *  dt  a  dt  u‘  ’ 


Equation  (2.17)  provides  a  means  of  computing  the  rate  of  change  of  cr„  from  3 h/bL 
as  long  as  we  can  estimate  all  the  other  terms  on  the  rhs.  Experimentation  has 
shown  that  the  last  term  is  always  at  least  two  orders  of  magnitude  smaller  than 
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the  first  two  terms  on  the  rhs,  so  that  term  is  dropped.  The  remaining  information 
is  readily  available  from  the  model  since  the  boundary-layer  temperature  can  be 
used  to  find  u  and  3ir/3f  can  be  computed  using  (2.7).  The  tendency  produced  by 
(2.17)  is  used  on  the  current  timestep  where  required  [for  example  in  (2.8)],  and  it  is 
also  used  to  calculate  the  new  value  of  cr ^  in  a  leapfrog  timestepping  scheme  (with 
Asselin  filter)  identical  to  that  used  for  the  other  prognostic  va-iables. 

e.  Boundary-Layer  treatment  of  grid  boxes  with  both  land  and  water 

The  presence  of  both  land  and  water  in  a  grid  element  requires  special 
consideration  in  the  boundary-layer  parameterization.  If  grid  boxes  were  treated  as 
either  all  water  or  all  land,  the  parameterization  could  be  coded  with  a  simple 
branch  to  handle  each  condition.  We  choose  to  allow  fractional  land-water  coverage 
in  grid  boxes  to  more  accurately  represent  the  actual  land  form  (see  section  4). 

We  assume  that  the  temperature  of  the  water  does  not  change  with  time 
during  the  integration  of  the  model.  Since  the  model  is  designed  to  run  for  only  24 
hours  after  initialization,  this  assumption  is  reasonable.  The  ten.perature  of  the 
ground  changes  quite  rapidly  durinp  a  given  integration,  and  this  temperature  is 
stored  in  the  variable  TG.  We  define  an  ‘‘effective  ground  temperature”  (ETG)  for 
each  grid  box  for  use  in  flux  calculations.  For  “all  lend”  grid  boxes,  the  effective 
ground  temperature  is  simply  TG.  For  “all  water”  grid  boxes,  ETG  is  the  water 
temperature,  TWATER.  For  partial  land  boxes,  the  effective  ground  temperature 
for  use  in  the  rest  of  the  parameterization  is  taken  to  be  a  weighted  average  of  TG 
and  TWATER  based  on  the  percent  land  coverage  in  the  grid  box.  The  quantity  TG 
is  allowed  to  change  as  usual,  since  land  surfaces  would  be  expected  to  heat 
normally  despite  the  presence  of  water  in  the  same  grid  box. 

Both  latent  and  sensible  heat  fluxes  are  affected  directly  by  the  presence  of 
a  water  surface.  The  sensible  heat  flux  is  driven  by  a  temperature  difference 
between  the  ground  or  wat°r  surface  and  the  air  just  above.  If  this  difference  is 
minimized  by  the  presence  of  constant-temperature  water,  the  sensible  heat  flux  will 
be  reduced.  The  presence  of  the  water  surface  also  ensures  an  endless  supply  of 
moisture  whicn  can  be  evaporated.  Thus,  fhe  latent  heat  flux  over  a  partial-land 
grid  box  will  tend  to  stay  higher,  and  this  will  prevent  the  land  from  heating  up  as 
fast  as  it  otherwise  would.  In  a  clear-sky,  daytime  situation,  the  net  effect  of 
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having  a  grid  box  partially  or  wholly  covered  by  water  is  to  reduce  the  growth  of 
the  convective  boundary  layer  and  slow  the  heating  of  the  effective  ground 
surface.  Figure  2  shows  the  height  of  the  convective  boundary  layer  for  a  range 
of  land/water  mixes.  The  solid  curve  represents  the  results  for  an  all-land  grid  box 
in  which  the  boundary  layer  grows  rapidly  after  sunrise.  Moderate  water  coverage 
in  the  grid  box  results  in  a  significant  reduction  of  the  boundary-layer  growth,  and 
for  the  0%  land  surface  case,  the  boundary  layer  grows  only  slightly.  For  the 


BOUNDARY  LAYER  HEIGHT 


6  12  18 
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Fig.  2.  Boundary  layer  height  as  a  function  of  time  for 
different  percentages  of  land  and  water  within  a  grid  box. 


—  1  im¬ 


partial  land  cases,  the  convective  boundary  layer  actually  collapses  into  a  stable 
one  between  1700  and  1800  local  time  (LST),  which  is  why  no  data  is  shown  at  18 
LST  for  these  runs. 

/.  Boundary  layer  transitions 

Transitions  between  stable  and  unstable  boundary  layer  regimes  have  been 
incorporated  into  the  boundary-layer  parameterization.  In  the  real  atmosphere, 
these  transitions  are  almost  certainly  abrupt  and  characterized  by  poorly  defined 
structures.  Within  the  model,  such  a  situation  would  lead  to  numerical  problems. 
In  addition,  there  is  no  known  parameterization  for  a  boundary  layer  in  transition, 
so  we  have  to  provide  a  mechanism  for  an  orderly  transition. 

The  transition  stage  is  determined  by  the  sign  of  the  sensible  heat  flux. 
When  the  direction  of  the  sensible  heat  flux  is  incompatible  with  the  nature  of  the 
boundary  layer,  the  transition  stage  is  set.  The  boundary  layer  height  is 
constrained  to  fall  from  its  current  location  to  a  height  of  90  m  (arbitrarily  chosen) 
in  15  min.  During  this  transition,  the  potential  temperature  structure  in  the 
boundary  layer  is  changed  as  little  as  possible,  thereby  allowing  the  proper 
parameterization  to  begin  with  the  current  structure  as  though  the  transition  had 
been  instantaneous.  The  potential  temperature  at  the  top  of  the  boundary  layer  is 
linearly  interpolated  between  its  current  value  and  the  potential  temperature  of  the 
ground.  Boundary  layer  moisture  is  conserved  during  this  process  of  boundary 
layer  height  fall,  again  to  preserve  current  conditions  as  much  as  possible.  When 
the  boundary  layer  height  reaches  SO  m  the  transition  ends.  The  boundary  layer 
then  begins  to  grow  again  according  to  the  stability  at  that  time. 

Transitions  from  unstable  to  stable  and  vice  versa  have  been  tested  and 
found  to  be  relatively  smooth  and  reliable.  Although  most  of  these  transitions  will 
take  place  at  dawn  or  dusk,  it  is  conceivable  that  a  major  change  in  cloud  cover  or 
air  mass  might  initiate  such  a  change  at  any  time  during  an  integration.  The  only 
aspect  of  this  which  will  not  be  handled  well  by  the  parameterization  is  a  situation 
in  which  a  well-mixed  convective  boundary  layer  grows  to  a  substantial  depth,  then 
the  radiation  is  interrupted  by  clouds,  causing  a  transition  to  a  stable  boundary 
layer,  followed  by  a  breakup  of  the  clouds  and  e  return  to  unstable  stratification. 
The  new  convective  boundary  layer  will  start  over  and  the  atmosphere  above  the 
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boundary  layer,  which  in  the  real  atmosphere  would  be  almost  dry  adiabatic  in 
structure,  will  be  slightly  stable  in  the  model.  This  is,  of  course,  a  problem  with 
the  limited  vertical  resolution  of  the  model,  not  the  fault  of  the  parameterization. 
Only  extensive  testing  will  be  able  to  show  if  this  is  a  serious  problem.  We  suspect 
that  most  cloud  cover  will  only  reduce  the  sensible  heat  flux  in  magnitude  rather 
than  actually  reversing  the  sign.  Hence,  the  boundary  layer  will  grow  much  more 
slowly  but  will  not  go  through  a  transition. 

An  example  of  the  results  after  two  transitions  is  shown  in  Fig.  3.  The 
ground  temperature  for  a  run  using  a  one-dimensional  version  of  the  boundary  layer 
parameterizations  is  shown  after  30  hours.  Clearly,  the  transitions  are  smooth  for 
this  variable.  The  boundary  layer  becomes  stable  near  1800  LST,  remains  stable 
until  0800  LST  on  the  second  day,  and  is  unstable  until  the  end  of  the  simulation  on 

GROUND  TEMPERATURE 
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4°  _ 
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Fig.  3.  Ground  temperature  for  a  30  h  simulation  using  the  one- 
dimensional  version  of  the  boundary-layer  parameterizations. 
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the  second  day  (1200  LST).  Dawn  occurs  just  before  0500  LST,  and  the  ground 
immediately  begins  to  heat  up  at  this  point.  The  boundary  layer  remains  stable 
until  the  ground  heating  is  sufficient  to  bring  the  ground  temperature  above  the 
surface  air  temperature.  This  happens  around  0800  LST  on  the  second  day.  One 
difference  between  the  two  morning  patterns  is  that  the  ground  temperature  reaches 
a  much  higher  value  on  the  second  day.  This  occurs  because  the  ground  surface 
dries  out  during  the  first  day,  allowing  the  available  incoming  energy  on  the  second 
day  to  be  used  much  more  for  heating  the  ground  and  less  for  evaporating  soil 
moisture. 


3.  Simulations  with  the  ETAS  Model 

The  five-layer  ^-coordinate  model  described  above  is  referred  to  as  the 
ETA5  model.  This  section  will  present  a  series  of  test  simulations  designed  to 
examine  the  capabilities  of  the  model  formulation  and  verify  the  correctness  of  the 
coding.  For  all  the  simulations  presented  here,  ETAS  is  being  run  without  the 
boundary  layer  parameterization  package  and  with  no  moisture.  This  allows  the 
hydrodynamic  portion  of  the  model  to  be  tested  without  the  complexities  introduced 
by  the  physical  parameterizations.  The  thermodynamic  structure  taken  for  all 
simulations  shown  here  is  the  U.S.  Standard  Atmosphere  profile. 

For  all  simulations  discussed  in  sections  3a  and  3b,  only  the  momentum 
equations  contain  an  eddy  diffusion  tei  m  with  the  coefficient  K  =■  2  V  10a  m:s_i, 
while  K  —  0  for  all  thermodynamic  variables.  Further,  no  additional  diffusion  was 
set  near  the  CGM  and  FGM  interface  region.  Additional  smoothing  techniques  are 
discussed  in  section  3c.  The  ability  of  the  model  to  correctly  allow  a  variable 
height  boundary  layer  is  tested  in  section  3b  by  forcing  the  boundary  layer  height 
change  explicitly. 

a.  Mountain  lee-wave  experiments  with  a  constant  boundary-layer  height 

The  simulation  of  mountain  lee-waves  is  generally  considered  an  important 
test  of  a  mesoscale  model  (Anthes  and  Warner  1978;  Nickerson  et  al.  1986).  SC 
showed  that  the  major  features  of  lee-waves  could  be  simulated  even  with  the 
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coarse  vertical  resolution  used  in  the  present  model.  This  section  will  present  the 
results  of  several  simulations  to  confirm  this  and  to  demonstrate  the  structure 
taken  on  in  the  ETA5  model  for  boundary  layers  of  different  depths.  Following  SC, 
we  will  conduct  these  tests  with  the  full  three-dimensional  model  rather  than  a  two- 
dimensional  analog.  In  all  the  simulations  shown,  a  1  km  ridge  with  a  Gaussian  cross 
section  is  aligned  north-south  in  the  FGM  domain.  The  ridge  terminates  smoothly 
within  the  FGM  domain  and  does  not  extend  into  the  CGM  domain  (see  FGM 
horizontal  plots  showing  height  contours  below).  For  all  simulations  shown,  a 
20  m  s~‘  u-component  wind  is  initially  specified  at  all  levels  and  held  constant  on 
the  outer  CGM  lateral  boundaries.  It  should  be  noted  that,  although  not  presented 
below,  simulations  with  a  negative  u  velocity  and  simulations  with  the  mountain 
oriented  east-west  have  been  carried  out  to  isolate  potential  coding  errors. 

Figure  4  shows  a  vertical  cross  section  through  the  center  of  the  FGM 
domain  for  a  lee-wave  simulation  in  which  the  Coriolis  parameter  was  set  to  zero. 
Figure  4a  shows  the  potential  temperature  field  at  12  h  simulated  time,  while  Fig.  4b 
shows  the  u-component  velocities.  The  dots  plotted  in  the  figures  show  the 
locations  of  thermodynamic  grid  points  on  TJ  surfaces  in  the  model.  In  Fig.  4a,  a 
thm  dashed  line  shows  the  top  of  the  boundary  layer  (that  is,  the  <jh  surface). 
The  isopleths  in  Fig.  4b  are  subjectively  analyzed  because  of  the  limited  number  of 
layers. 

The  structure  shown  in  Fig.  4  matches  the  qualitative  structure  of  a 
mountain  lee-wave,  with  the  wave  structure  sloping  into  the  wind  with  height  and  a 
reversal  of  the  velocity  perturbation  directly  above  the  mountain  aloft  compared  to 
the  layer  near  the  surface.  The  boundary  layer  height  for  this  simulation  was  set 
by  specifying  crn  =  0.90.  This  results  in  a  boundary  layer  that  is  approximately  91 
mb  (about  790  m)  deep.  The  ^-coordinate  formulation  then  spaces  the  remaining 
four  layers  above  the  boundary  layer  evenly  in  ct.  Figure  5  shows  a  horizontal 
plot  for  this  same  simulation  and  time  showing  the  sea-level  pressure  (reduced  using 
a  constant  lapse  rate  hydrostatic  formula)  and  boundary  layer  winds.  The  weak 
upstream  damming  and  lee  trough  are  evident  in  the  pressure  contours  representing 
vaiues  above  and  below  the  initially  specified  sea-level  pressure  of  1013.25  mb. 
The  deflection  of  boundary  layer  flow  around  the  mountain  is  also  evident  in  the 
-vino  field.  The  time  required  by  the  model  to  reach  a  steady  flow  can  be  assessed 
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Fig.  4.  Vertical  cross  section  throught  the  FGM  domain  at  1*.  h  simulated  time  for 
lee  wave  test  with  no  Coriolis  effect  showing  (a)  the  potential  temperature,  and  (b) 
the  u-component  velocity.  Isentropes  are  labeled  in  K  along  the  right  in  (a)  and  grid 
point  values  of  the  velocity  are  plotted  in  (b)  in  m  s~’.  The  top  of  the  boundary 
layer  is  shown  in  (a)  by  a  thin  dashed  line  =  0.90). 


Fig.  5.  Horizontal  plot  for  simulation  shown  tn  Fig.  4.  Solid  contours  give 
sea-level  pressure  (in  mb),  thin  dashed  contours  are  for  surface  topography 
(in  in)  of  the  mountain  ridge,  and  boundary  layer  winds  (one  full  barb 
represents  10  m  s~l)  are  plotted  at  each  grid  point. 
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in  the  FGM  domain  for  the  simulation.  This  is  shown  in  Fig.  6.  After  a  period  of 
oscillations  of  rapidly  decreasing  amplitude  for  the  first  few  hours,  the  model 
simulation  becomes  quite  steady.  Despite  the  steadiness  of  the  solution,  some  small 
amplitude  2Ax  noise  develops  in  the  temperature  field  during  the  last  6  h  of  the 
simulation.  This  is  noticable  in  Fig.  4a  near  the  western  boundary  (the  left  side  of 
the  plot).  It  is  not  surprising  that  some  noise  would  develop  given  the  complete 
lack  of  smoothing  of  temperatures  in  the  model.  In  fact,  it  is  more  surprising  that 
the  model  could  be  integrated  for  12  h  with  a  fairly  large  amplitude  disturbance  and 
generate  so  little  noise. 

We  carried  out  an  experiment  similar  to  that  described  above  except  with  the 
Coriolis  effect  included  by  setting  /  =  lO-"  s-1  and  specifying  an  initial  pressure 
gradient  in  balance  with  a  20  m  s_1  wind.  Figures  showing  the  same  fields  as 
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Fig.  6.  Average  rate  of  change  of  ir  versus  time  for  FGM  domain  during 
the  simulation  shown  in  Figs.  4  and  5. 
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Fig.  8.  As  in  Fig.  5,  except  that  the  Coriolis  effect  was  included. 
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Figs.  4—6  are  shown  in  Figs.  7—9.  The  cross  section  plots,  Ftgs.  7a  and  7b,  are 
nearly  identical  to  those  produced  with  no  Coriolis  effect,  except  that  the  wave  has 
a  slightly  smaller  amplitude.  The  horizontal  plot,  however,  looks  considerably 
different  since  the  large  scale  pressure  field  is  now  in  near-geostrophic  balance  with 
the  winds.  A  southward  turning  of  the  wtnds  is  evident  due  to  a  slight  imbalance  of 
the  initially  specified  pressure  field  and  the  wind  field  after  its  rapid  adjustment  to 
the  mountain  barrier.  The  small  amplitude  noise  that  develops  late  in  the 
simulation  is  apparent  in  the  isobars  in  Fig.  8.  Despite  this  low  level  of  noise.  Fig. 
9  indic°*es  that  the  solution  is  quite  steady  after  the  initial  adjustment  period,  and 
this  is  confirmed  in  Fig.  10  which  shows  the  superimposed  sea-level  pressure  fields 
for  8  h,  10  h,  and  12  h.  For  completeness.  Fig.  11  shows  a  cross  section  of  potential 
temperature  for  the  CGM  domain.  The  solution  shown  in  the  central  region 


Fig.  9.  Average  rate  of  change  of  ir  versus  time  for  FGM  domain  during 
the  simulation  shown  in  Figs.  7  and  8. 
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Fig.  10.  Sea-level  pressure  fields  at  8  h,  10  h,  and  12  h,  for  the  simulation 
shown  in  Figs.  7-9.  Solid  contours  give  sea-level  pressure  (in  mb),  thin 
dashed  contours  are  for  surface  topography  (in  m)  of  the  mountain  ridge. 
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(between  the  vertical  dashed  lines)  is  actually  the  FGM  solution  averaged  to  CGM 
grid  points.  Note  the  smoothness  of  this  solution,  even  at  the  FGM/CGM  interface. 

Figure  12  shows  a  potential  temperature  cross  section  for  a  12  hour 
simulation  in  which  the  boundary  layer  was  specified  to  be  fairly  thin.  For  this 
simulation  a\  =  0.96,  which  leads  to  a  boundary  layer  thickness  of  about  36  mb  (or 
about  308  m).  In  all  other  respects  this  simulation  was  identical  to  that  shown  in 
Fig.  7.  The  change  in  boundary  layer  depth  affects  the  locations  of  all  levels  in 
the  model,  so  the  vertical  resolution  aloft  is  even  coarser  than  before.  The 
qualitative  features  of  the  solution  are  quite  similar  to  the  other  simulations,  but 
some  differences  are  worth  noting.  Keeping  the  lowest  layer  much  thinner  increases 
the  influence  of  the  mountain  on  the  boundary  layer  wind  field,  as  can  be  seen  by 
comparing  Fig.  13  with  Fig.  8.  For  the  thin  boundary  layer  case  (Fig.  13),  the 
upwind  deceleration  of  the  winds  is  more  pronounced,  resulting  in  a  sharper 
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Fig.  11.  Vertical  cross  section  of  potential  temperature  through  CGM 
domain  at  12  h  for  simulation  shown  in  Figs.  7-9.  The  vertical  dashed 
lines  represent  the  FGM/CGM  interface  location. 
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upwind  ridge,  and  there  is  more  lateral  deflection  of  the  winds  around  the  mountain 
ridge. 

Figures  14  and  IS  show  the  opposite  situation,  where  the  boundary  layer 
height  was  set  to  be  thicker  than  in  previous  runs.  Here,  <7h  =  0.80,  which  results 
in  a  boundary  layer  thickness  of  about  181  mb  (about  1645  m).  For  this  value  of 
the  model’s  five  layers  are  equally  spaced  in  a.  The  discussion  of  nie  previous 
paragraph  holds  in  reverse  now,  as  the  mountain  has  less  influence  on  the  boundary 
layer  winds  since  the  boundary  layer  top  is  now  well  above  the  top  of  the  ridge.  It 
is  interesting  to  note  that  the  2Ax  noise  is  somewhat  greater  for  the  thick 
boundary  layer  than  it  is  for  the  thin  one. 

Since  there  is  no  imposed  (or  implied)  capping  inversion  of  the  boundary 
layers  in  these  simulations,  the  differences  in  the  solutions  for  the  different  crh 
values  are  purely  an  artifact  of  the  lack  of  sufficient  vertical  resolution.  Given 


Fig.  12.  Vertical  cross  section  of  potential  temperature  through  FGM 
domain  at  12  h  showing  potential  temperature  for  thin  boundary  layer 
simulation  (cth  =  0.96). 
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Fig.  13.  Horizontal  plot  of  FGM  domain  at  12  h  for  thin  boundary  layer 
simulation  (cr„  =  0.96).  Solid  contours  give  sea-level  pressure  (in  mb),  thin 
dashed  contours  are  for  surface  topography  (in  m)  of  the  mountain  ridge, 
and  boundary  layer  winds  (one  full  barb  represents  10  m  s_l)  are  plotted  at 
each  grid  point. 
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many  layers,  the  model  would  resolve  more  flow  around  the  mountain  in  the  lowest 
layers  and  somewhat  less  in  the  layers  just  above  the  surface.  The  current  model 
must  be  viewed  as  providing  a  sense  of  the  average  flow  over  the  layer  defined 
here  as  the  boundary  layer.  Since  each  solution  represents  a  different 
approximation  of  this  type,  each  solution  will  present  slightly  different  fields  of 
winds,  pressure  and  temperatures.  The  fact  that  the  solutions  are  very  similar 
over  a  wide  range  of  boundary  layer  depths  lends  support  to  the  ability  of  the  77- 
coordinate  system  to  represent  the  real  atmosphere. 


Fig.  14.  Vertical  cross  section  through  FGM  domain  at  12  h  showing 
potential  temperature  for  thick  boundary  layer  simulation  (ct„  «■*  0.80). 
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Fig.  15.  Horizontal  plot  of  FGM  domain  at  12  h  for  thick  boundary  layer 
simulation  (cr.,  =  0.80).  Solid  contours  give  sea-level  pressure  (in  mb),  thin 
dashed  contours  are  for  surface  topography  (in  m)  of  the  mountain  ridge, 
and  boundary  layer  winds  (one  full  barb  represents  10  m  s'1)  are  plotted  at 
each  grid  point. 


28  — 


b.  Experiments  with  a  changing  boundary-layer  height 

The  simulations  presented  in  the  previous  section  all  treated  the  boundary 
layer  height,  ah,  as  a  constant  for  the  entire  period  of  the  integration.  Inspection 
of  the  equation  set  shows  that  if  dcr^/dt  =  0,  the  model  reduces  to  a  cr-coordinate 
model.  It  is  therefore  important  to  test  the  ability  of  the  model  to  maintain  the 
proper  atmospheric  structure  even  when  the  boundary  layer  height  is  changing.  A 
simple  example  of  this  is  shown  in  Figs.  16a  and  16b.  These  show  the  results  of 
two  separate  integrations  after  4  h  simulated  time,  with  flat  terrain  and  a  u- 
component  velocity  of  20  m  s_1  at  all  levels.  Figure  16a  has  a  constant  boundary 
layer  height  at  <y„  =  0.96,  and  shows  essentially  no  change  during  the  integration. 
For  Fig.  16b,  however,  the  model  was  initiallized  with  ah  —  0.96  throughout  the 
domain,  but  the  boundary  layer  height  was  artificially  forced  upward  in  the  central 
region  of  the  FGM  by  specifying  dh/dt  as  a  Gaussian  forcing  function  whose 
central  value  was  400  m  h~l.  As  can  be  seen  in  Fig.  16b,  this  resulted  in  the 
dynamic  adjustment  of  the  model  layers  to  the  changing  boundary  layer  thickness, 
but  the  isentropic  surfaces  remained  nearly  horizontal  —  as  they  should.  An 
interesting  aspect  of  the  ^-coordinate  model  is  that  even  though  the  flow  is 
horizontal  in  physical  space,  because  of  the  upward  bulge  of  the  T]  surfaces  there  is 
a  significant  flow  through  the  surfaces  that  represents  a  vertical  velocity  in  the  T] 
system.  Figure  17  shows  the  geostrophic  adjustment  for  both  simulations  in  terras 
of  the  rate  of  change  of  surface  pressure  (note  the  scales  are  expanded  compared  to 
previous  figures  of  this  type).  As  can  be  seen,  the  movement  of  the  crH  surface 
has  no  effect  on  the  adjustment  process  in  the  model. 

A  somewhat  different  experiment  was  carried  out  by  forcing  a  rising 

boundary  layer  during  a  lee-wave  simulation.  In  this  case,  the  simulation  started 
with  ah  =  0.96  throughout  the  domain  and  these  conditions  were  held  for  6  h 
simulated  time.  Then,  the  boundary  layer  was  forced  to  rise  by  specifying  dh/dt 
*=  400  m  s-1  at  all  grid  points.  This  rate  of  rise  was  held  until  =  0.80  when 
the  rate  of  rise  was  set  back  to  zero  (at  about  9.3  h).  Figure  18  shows  a  cross 
section  for  this  simulation  at  12  h.  Comparison  with  Fig.  14  shows  that  the  growth 
of  the  boundary  layer  did  little  to  affect  the  final  solution,  despite  the  large 
change  in  physical  elevation  of  the  model  surfaces  and  the  implied  vertical 

velocities  that  went  along  with  the  rising  surfaces.  Because  the  flow  relaxation 

boundary  condition  applied  on  the  lateral  boundaries  of  the  CGM  holds  the 
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boundary  values  fixed,  these  boundary  points  retained  a  boundary  layer  height  of 
cr^  =  0.96  throughout  the  simulation  —  requiring  an  additional  adjustment  process 
to  take  place  in  the  CGM  collar  surrounding  the  FGM  domain.  A  cross  section 
through  the  CGM  which  shows  this  is  shown  in  Fig.  19.  There  is  some  noise  in  and 
just  above  the  boundary  layer  inside  the  relaxation  region  near  the  outer  lateral 
boundaries,  but  the  solution  is  really  quite  well-behaved.  Figures  20,  21,  and  22 
show  horizontal  plots  of  the  FGM  sea-level  pressure  field  at  6  h,  8  h,  and  12  h, 
respectively.  Figures  20  and  22  look  almost  identical  to  Figs.  13  and  15,  as  they 
should  since  they  represent  lee-wave  solutions  for  the  same  vertical  structure  of 
the  model.  Figure  21  shows  the  solution  while  the  boundary  layer  is  growing,  and 
although  it  is  a  little  noisier,  the  solution  is  again  quite  well-behaved.  In  addition 
to  the  plots  shown  above,  the  lack  of  impact  on  the  model’s  ability  to  simulate 


TIME  (hr) 

Fig.  17.  Average  rate  of  change  of  it  versus  time  for  FGM  domain  during 
the  simulations  shown  in  Figs.  16a  (solid)  and  16b  (dashed).  Note  that 
dashed  line  only  deviates  from  solid  line  at  the  end  of  the  simulation. 
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Fig.  18.  Vertical  cross  section  through  FGM  domain  at  12  h  showing 
potential  temperature  for  rising  boundary  layer  simulation. 


Fig.  19.  Vertical  cross  section  through  CGM  domain  at  12  h  showing 
potential  temperature  for  rising  boundary  layer  simulation. 
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Fig.  20.  Horizontal  plot  of  FGM  domain  at  6  h  for  growing  boundary  layer 
simulation  (ct*  =  0.96  at  this  time).  Solid  contours  give  sea-level  pressure 
(in  mb),  thin  dashed  contours  are  for  surface  topography  (in  m)  of  the 
mountain  ridge. 
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Fig.  21.  Horizontal  plot  of  FGM  domain  at  8  h  for  growing  boundary  layer 
simulation  ( a H  =  0.86  at  this  time  and  is  rising  at  400  m  s_l).  Solid 
contours  give  sea-level  pressure  (in  mb),  thin  dashed  contours  are  for 
surface  topography  (in  mi. 
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Fig.  22.  Horizontal  plot  of  FGM  domain  at  12  h  for  growing  boundary 
layer  simulation  (crh  =  0.80  at  this  time).  Solid  contours  give  sea-level 
pressure  (in  mb),  thin  dashed  contours  are  for  surface  topography  (in  m). 
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the  lee-wave  phenomenon  despite  a  drastic  change  in  boundary  layer  height  is  shown 
in  Fig.  23.  Here  the  average  rate  of  change  of  ir  is  plotted  for  the  simulation 
discussed  here.  While  a  very  small  adjustment  in  the  model  is  apparent  between  6  h 

and  12  h  while  the  boundary  layer  is  rising,  it  is  clear  that  the  model  equations 

allow  the  boundary  layer  change  to  occur  without  disturbing  the  geostrophic 
balance  within  the  model. 

Several  other  simulations,  which  are  not  shown  here,  have  been  carried  out 
in  which  the  boundary  layer  was  artificially  raised  or  lowered  to  test  the 
robustness  of  the  77-formulation.  In  all  cases,  the  model  performed  well.  This 
includes  an  extreme  case  in  which  the  Gaussian  forcing  function  was  imposed  on  a 
mountain  lee-wave  simulation  for  the  full  12  h  so  that  the  boundary  layer  directly 
over  the  ridge  grew  to  almost  9  km  in  depth.  Even  in  this  case  the  model  solution 

remained  stable,  although  the  solution  was  not  a  vary  good  simulation  of  lee-waves 


Fig.  23.  Average  rate  of  change  of  ir  versus  time  for  FGM  domain  during 
the  growing  boundary  layer  simulation  shown  in  Figs.  18—22. 
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since  the  lee-wave  could  not  be  resolved  by  the  single  very  thick  layer  just  above 
the  mountain. 

c.  Eddy  diffusion  sensitivity  tests 

All  of  the  simulations  shown  in  this  section  were  carried  out  with  no 
diffusion  on  temperature  and  no  additional  smoothing  operators  on  any  variables 
near  the  FGM/CGM  interface.  This  was  desirable  in  early  testing  so  that  even 
small  incompatibilities  between  CGM  and  FGM  solutions,  errors  in  coding,  or 
accumulating  errors  due  to  a  nonconstant  boundary  layer  depth  would  be  readily 
apparent.  It  is  clear  from  inspection  of  Fig.  8,  for  example,  that  some  noise  is 
present  in  the  model  that  should  be  removed  by  smoothing.  Figure  10  shows  that 
there  are  actually  two  types  of  noise  present.  On  the  eastern  boundary  is  a  steady 
2&x  noise  that  is  related  to  the  overspecification  of  pressure  on  the  outflow 
FGM/CGM  interface  and  its  impact  on  the  velocities.  We  will  call  this  the  “outflow 
noise.”  The  western  half  of  the  domain,  upwind  of  the  mountain  ridge  also  has  some 
2Az  noise  which  appears  more  transient.  This  appears  to  be  simply  the  noise  that 
is  typical  in  primitive  equation  models  when  the  nonlinear  advective  terms  are 
present,  and  is  related  to  the  lack  of  any  smoothing  on  temperature.  This  noise  is 
easily  removed  by  a,  ding  even  a  very  small  amount  of  diffusion  to  the  temperature. 
Figure  24  shows  a  horizontal  plot  of  a  lee-wave  simulation  with  uh  =  0.80  in  which 
an  eddy  diffusion  coefficient  of  K  =  1  X104  nFs~‘  was  specified  for  temperature. 
This  value  is  only  0.05  times  that  used  for  the  momentum  terms.  We  chose  to  use 
Ok  =  0.80  for  this  simulation  because  the  noise  is  somewhat  more  apparent  in  the 
thick  boundary  layer  simulations  (see  Figs.  14  and  15).  Note  that  while  the  extra 
diffusion  on  temperature  removed  the  regular  noise  in  the  western  portion  of  the 
domain,  it  did  not  affect  the  “outflow  noise"  on  the  eastern  boundary. 

Zhang  et  al.  (1986)  used  both  the  Newtonian  damping  scheme  discussed  in 

section  2  and  a  larger  eddy  diffusion  coefficient  near  the  FGM/CGM  interface  to 

control  the  noise  produced  by  the  overspecification  of  pressure  on  outflow 

boundaries.  Figure  25  shows  a  simulation  (again  with  crh  =  0.80)  in  which  there  is 
no  diffusion  on  temperature  but  the  Newtonian  damping  scheme  is  used  on  the 

momentum  terms  (with  Td  -  20AG.  While  this  reduces  the  outflow  noise,  it  does 
not  eliminate  it.  Additional  tests  (not  shown)  indicated  that  additional  eddy 
diffusion  was  required  near  the  boundary  to  remove  this  noise.  As  discussed  in 
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Fig.  24.  Hor.zontal  plot  of  FGM  domain  at  12  h  for  simulation  with  weak 
diffusion  on  temperature  (crH  =  0.80!.  Solid  contours  give  sea-level 
pressure  (in  mb),  thin  dashed  contours  are  for  surface  topography  (in  m)  of 
the  mountain  ridge,  and  boundary  layer  winds  (one  full  barb  represents 
10  m  s"‘)  are  plotted  at  each  grid  point. 
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Fig.  25.  Horizontal  plot  of  FGM  domain  at  12  h  for  simulation  with 
Newtonian  damping  on  momentum  terms  at  the  CGM/FGM  interface 
(crn  =  0.80).  Solid  contours  give  sea-level  pressure  (in  mb),  thin  dashed 
contours  are  for  surface  topography  (in  m)  of  the  mountain  ridge,  and 
boundary  layer  winds  (one  full  barb  represents  10  m  s-1)  are  plotted  at 
each  grid  point. 


39  — 


section  2,  the  Newtonian  damping  scheme,  as  applied  here,  reduces  to  a  simple  filter 
and  is  nearly  identical  to  increased  diffusion.  We  therefore  decided  to  simply 
control  this  noise  by  increasing  the  diffusion  near  the  boundary. 

The  combination  of  these  two  noise  control  techniques  results  in  simulations 
such  as  that  shown  in  Fig.  26.  For  this  simulation,  the  eddy  diffusion  on 
temperature  is  included  (as  in  Fig.  24),  and  the  diffusion  on  all  terms  is  modified 
near  the  FGM/CGM  interface  to  be  twice  the  interior  value  2  grid  intervals  away 
from  the  interface  and  three  times  the  interior  value  1  grid  interval  from  the 
interface.  This  modest  additional  smoothing  affects  the  simulation  in  only  minor 
ways  and  greatly  reduces  the  noise.  Some  of  the  modification  of  the  isobar  pattern 
near  the  mountain  is  a  result  of  the  small  accumulating  error  discussed  by  SC  which 
occurs  whenever  diffusion  is  applied  to  temperature  on  a  or  r|  surfaces  in  complex 
terrain.  The  very  small  value  of  the  diffusion  parameter  makes  this  error  small 
enough  to  provide  acceptable  results  even  in  this  very  steady  flow,  and  it  should 
be  unnoticeable  in  a  more  transient  flow. 

Figures  27  and  28  provide  the  results  for  a  simulation  of  the  model  with  the 
added  smoothing  for  ah  —  0.90.  These  may  be  compared  directly  with  Figs.  7  and 
8.  While  there  are  some  subtle  differences,  the  simulation  with  the  extra  diffusion 
is  a  qualitatively  good  solution  for  mountain  lee-waves  and  has  very  little  noise. 
We  feel  that  Figs.  27  and  28  represent  the  model  as  it  will  be  run  from  now  on. 


4.  Domain  for  Prototype  Testing 

An  important  aspect  of  this  study  is  the  testing  of  a  prototype  model  using 
real  data  in  clear  weather  and  stratiform  precipitation  situations.  Such  testing 
requires  a  specific  geographic  location.  Ideally,  this  location  should  experience  a 
wide  variety  of  weather  phenomena  such  as  sea-breeze  circulations,  synoptic  scale 
stratiform  precipitations  events,  orographic  enhancement  of  precipitation,  and  other 
mesoscale  patterns  for  which  the  model  is  intended  to  provide  guidance.  Southern 
New  England  represents  one  of  several  locations  that  could  be  considered,  and  given 
our  obvious  desire  for  operational  guidance  in  our  own  area,  we  have  chosen  this 
geographic  region  for  testing  of  the  prototype  model. 
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Fig.  26.  Horizontal  plot  of  FGM  domain  at  12  h  for  simulation  with  weak 
diffusion  on  temper  lture  and  extra  diffusion  near  FGM/CGM  interface 
(crh  =*  0.80).  Solid  contours  give  sea-level  pressure  (in  mb),  thin  dashed 
contours  are  for  surface  topography  (in  m)  of  the  mountain  ridge,  and 
boundary  layer  winds  (one  full  barb  represents  10  m  s-1)  are  plotted  at 
each  grid  point. 
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Fig.  27.  As  in  Fig.  7,  except  for  the  ad< 
diffusion  near  the  boundaries. 
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Fig.  28.  As  in  Fig.  8,  except  for  the  added  weak  diffusion  on  temperature 
and  extra  diffusion  near  the  boundaries. 
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The  domain  for  the  fine  grid  mesh  is  shown  in  Fig.  29  (only  the  locations  of 
“thermodynamic"  grid  points  are  shown).  The  coarse  grid  mesh  which  surrounds  the 
FGM  domain  is  shown  in  Fig.  30.  In  addition  to  being  a  coastal  location,  this  area 
has  complex  terrain  which  is  rich  in  detail  on  meao-0  scales.  This  is  quite  evident 
in  the  contour  plot  of  the  topography  shown  in  Fig.  31.  This  plot  was  constructed 
by  reading  the  heights  of  FGM  points  directly  from  U.S.  Geologic  Survey 
topographic  maps.  The  field  is  quite  noisy  at  this  20  km  resolution,  so  it  is  clear 
some  processing  is  required.  We  anticipate  that  some  form  of  “envelope"  orography 
will  provide  the  best  overall  representation  of  the  effects  of  topography  on  the 
flow  (Wallace  et  al.  1983).  In  order  to  compute  an  envelope  field,  the  terrain 
heights  were  recorded  for  nine  equally  spaced  points  within  each  grid  box.  A  full 
envelope  is  considered  to  be  the  mean  of  these  nine  points  plus  2.0  times  the 

standard  deviation  of  them  (we  will  refer  to  this  multiplicative  factor  as  the 

“envelope  parameter”).  A  full  envelope,  with  the  envelope  parameter  equal  to  2.0, 
results  in  some  grid  point  elevations  being  considerably  higher  than  any  of  the  raw 
grid  point  heights.  Therefore,  we  anticipate  using  a  “partial  envelope”  in  which  the 
envelope  parameter  is  set  to  1.0.  As  shown  in  Fig.  32,  this  partial  envelope 

provides  a  desirable  smoothing  of  the  terrain  without  overemphasizing  any 

particular  feature.  Several  mesoscale  orographic  features  thought  to  be  important 
to  weather  in  the  New  England  area  are  obvious  in  the  smoothed  fields  of  Fig.  32. 
These  include  the  Adirondack  mountains,  the  Hudson  River  Valley,  the  Green 
Mountains,  the  Berkshires,  the  Connecticut  River  Valley,  the  White  Mountains,  and 
the  Worcester  Plateau  (see  Fig.  33  for  a  map  showing  the  locations  of  these 
features).  A  potential  problem  is  that  Mt.  Washington  falls  almost  directly  on  a 
grid  point  at  the  northern  edge  of  the  FGM  domain  and  therefore  tends  to  dominate 
the  terrain  field  in  northern  New  Hampshire.  This  particular  point  may  need  to  be 
reduced  in  height  in  order  to  avoid  numerical  problems  near  the  FGM  boundary. 
Further  processing  of  the  terrain  data  includes  an  adjustment  to  the  fields  in  the 
region  where  the  CGM  and  FGM  overlap  so  that  coincident  FGM  and  CGM  grid 
points  have  the  same  elevation  as  described  by  Zhang,  et  al.  (1986).  Further,  an 
additional  smoother  will  be  applied  to  the  terrain  field  to  remove  any  large 
amplitude  2Ax  component  which  may  tend  to  excite  noise  in  the  model. 

As  discussed  in  section  2,  the  model’s  boundary  layer  parameterization 
scheme  allows  for  a  grid  box  to  be  composed  of  a  mixture  of  land  and  water.  This 
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Fig.  29.  FGM  domain  showing  location  of  “thermodynamic-’  grid  points. 
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Fig.  30.  CGM  domain  showing  location  of  “thermodynamic”  grid  points. 
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Fig.  31.  Unsmoothed  model  topography  for  FGM  domain  (contour  interval 
is  50  m). 
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Fig.  32.  FGM  topography  for  a  partial  envelope  (envelope  parameter  is  1.0, 
contour  interval  is  50  m). 
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Fig.  33.  Map  of  FGM  domain  showing  major  geographical  features 
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means  the  coastline  need  not  be  considered  a  "blocky”  approximation  to  the  actual 
coastline  shown  in  Fig.  29,  and  it  also  means  that  the  influence  of  inland  bodies  of 
water  can  be  included  in  the  model-derived  fluxes.  The  U.S.  Geologic  Survey  maps 
were  again  used  to  make  a  careful  subjective  estimate  of  the  percent  of  water 
coverage  for  each  grid  box,  centered  on  the  thermodynamic  grid  points  shown  in  Fig. 
29.  A  contour  plot  of  percent  water  coverage  for  the  FGM  domain  is  shown  in  Fig. 
34.  Notice  how  well  even  subtle  variations  in  the  coastline  are  captured  in  the 
contouring,  and  the  influence  of  some  important  inland  bodies  of  water,  such  as  the 
Hudson  River,  the  Quabbin  Reservoir,  and  Lake  Winmpesaukee  (see  Fig.  33).  The 
island  of  Nantucket  presented  a  problem  in  this  analysis.  Review  of  Fig.  29  shows 
that  Nantucket  lies  nearly  in  the  middle  of  a  square  formed  by  four  grid  points. 
This  means  that  its  land  is  almost  equally  divided  between  the  four  grid  boxes 
represented  by  those  four  grid  points.  Rather  than  have  the  model  perceive  a 
40  ,<40  km  area  of  partial  land  coverage,  the  decision  was  made  to  “move”  Nantucket 
approximately  10  km  northwest  so  that  its  land  mass  fell  entirely  within  one  grid 
box.  This  move  is  obvious  in  Fig.  34,  in  which  the  partial  land  contour  representing 
Nantucket  is  displaced  northwestward  of  the  geographic  location  of  the  island. 


5.  Comparison  of  the  boundary-layer  parameterizations  with  the  O’Neill  dataset 

A  series  of  comparison  runs  have  be  in  made  with  the  parameterizations  being 
run  in  a  one-dimensional  form.  For  this  report,  we  shall  show  comparisons  with  data 
taken  during  the  O’Neill  boundary  layer  experiment  in  Kansas  (Lettau  and  Davidson 
1957). 

Figure  35  shows  the  comparison  graphs  for  the  various  fluxes.  In  the  top 
panel,  the  net  radiation  (NR)  and  soil  heat  flux  (GS)  are  shown.  The  current  model 
results  are  represented  by  the  small  circles  connected  with  solid  lines.  The  squares 
connected  with  dashed  lines  are  model  results  from  Colby  (1983),  which  are  shown 
for  comparison  because  the  current  model  was  adapted  from  these  previous 
parameterizations.  The  data  from  the  O’Neill  experiment  appears  as  crosses  with  no 
lines  connecting  them.  The  results  for  NR  compare  very  well  with  both  the 
previous  modeling  results  and  the  observations.  The  current  model  underforecasts 
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Fig.  34.  Contour  plot  of  FGM  domain  showing  the  percentage  of  water  in 
each  grid  box  (contour  interval  is  15%) 
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NTT  RADIATION  SNR)  AND 
SOIL  HEAT  FLUX  (GS) 
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Fig.  35.  Various  fluxes  for  O’Neill  case  data  from  the  one-din  ensional  version  of 
the  boundary-layer  parameterization  (solid  line  connecting  dots),  previous  model 
(dashed  line  connecting  squares),  and  observations  (crosses). 
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the  peak  NR  by  about  5%,  compared  to  an  apparent  overforecast  of  about  10%  by 
the  previous  model  formulation. 

The  sensible  and  latent  heat  fluxes  (SH,  LH)  are  also  shown  in  Fig.  35.  It  is 
obvious  that  the  flux  parameterization  is  successful  in  reproducing  the  observed 
values  within  10%  of  the  actual  measurements.  Although  there  appears  to  be  a 
large  difference  between  model  and  observations  for  the  latent  heat  flux,  the 
observations  are  subject  to  large  measurement  errors.  The  magnitude  of  these 
errors  is  ±  50  meal  min-1;  clearly  both  models  reproduce  the  observations  within 
the  expected  error. 

The  soil  heat  flux  (GS)  is  not  forecast  well  by  either  model  in  the  early 
morning,  with  the  current  model  overforecasting  more  than  the  previous  model. 
Both  models  rapidly  converge  with  the  measurements,  however,  between  0900  and 
1000  LST.  After  this,  both  models  coincide  with  the  measurements.  The  early 
morning  overforecast  of  GS  is  probably  responsible  for  the  subsequent  forecast 
errors  in  ground  temperature  (TG)  evident  in  Fig.  36.  It  is  apparent  in  this  figure 
that  the  ground  temperature  is  predicted  to  be  too  high  by  about  1  K  during  the 
early  morning.  Later,  this  error  reverses  sign  when  the  measured  ground 
temperature  exceeds  the  forecasts  of  both  models.  The  current  model  is  especially 
low,  by  about  3.5  K  near  1200  LST.  The  exact  reason  for  this  model  discrepency, 
particularly  the  difference  between  models,  is  not  known.  We  plan  to  investigate 
this  further. 

Figure  36  illustrates  the  rise  of  the  inversion  during  the  day,  showing  the 
height  of  the  top  of  the  boundary  layer  in  meters.  The  current  model  seems  to 
force  a  too-rapid  growth  of  the  boundary  layer,  being  at  times  almost  250  m  higher 
than  observed.  This  graph  clearly  demonstrates  the  limitation  imposed  by  the 
limited  vertical  resolution  of  the  model.  In  the  real  atmosphere,  as  well  as  the 
previous  model,  the  inital  boundary  layer  growth  was  suppressed  by  a  very  strong 
stable  layer.  Once  this  layer  was  eroded,  the  atmosphere  above  was  only  marginally 
stable,  allowing  rapid  growth  —  as  shown  by  the  previous  model  results  between 
1400  and  1700  LST.  In  the  current  model,  this  two-layer  stability  is  averaged  into 
one  stable  layer,  allowing  too-rapid  growth  in  the  morning  and  too-slow  growth  in 
the  afternoon.  After  1700  LST,  the  distinction  between  the  boundary  layer  and  the 
layer  immediately  above  disappeared.  In  the  real  atmosphere,  the  boundary  layer 
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Fig.  36.  As  in  Fig.  35,  except  for  height  of  inversion,  ground  temperature,  and 
“surface"  temperature. 
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grew  rapidly  into  the  layer  above,  which  is  what  both  the  previous  and  the  current 
models  demonstrate. 

We  believe  that  the  growth  problem  and  the  ground  temperature  problem 
come  together  in  the  forecast  of  the  surface  temperature  (TS)  shown  in  Fig.  36. 
The  error  between  the  current  model  and  the  observations  is  nearly  4  K  at  its 
worst.  We  plan  further  comparison  runs  to  determine  whether  this  magnitude  of 
error  will  be  a  common  occurrence. 


6.  Conclusions  and  future  work 

The  work  outlined  here  demonstrates  that  the  H-coordinate  model  we  have 
developed  is  robust  and  can  model  the  real  atmosphere  reasonably  well  despite 
having  a  very  limited  vertical  resolution.  In  addition,  the  boundary-layer 
parameterizations  work  quickly  and  well,  forecasting  the  gross  characteristics  of 
both  stable  and  unstable  boundary  layers  with  good  accuracy. 

The  next  tasks  which  will  be  undertaken  are  (1)  the  inclusion  of  the 
boundary-layer  parameterizations  into  the  ^-coordinate  model,  and  (2)  the  testing  of 
the  model  in  both  idealized  and  real-world  situations.  Inclusion  of  the  boundary- 
layer  code  is  somewhat  tricky  since  some  of  the  calculations  must  be  performed  in 
separated  subroutines  for  consistency.  Many  of  the  calculations  done  in  the  one¬ 
dimensional  form  of  the  boundary-layer  parameterizations  will  be  removed  (for 
example,  the  budget  calculation  used  to  find  the  boundary  layer  mixing  ratio),  and 
the  output  of  results  must  be  reworked  to  provide  the  correct  variables  to  the  main 
model. 

When  we  are  sure  that  the  boundary-layer  parameterizations  are  working 
correctly  in  the  full  model,  we  plan  to  run  several  tests  of  the  model.  One  of  the 
first  tests  will  be  a  simple  idealized  sea-breeze  circulation.  Other  tests  will  involve 
diurnal  variations  of  flow  over  simplified  terrain,  similar  to  the  ridge  used  in  the 
lee-wave  simulations.  These  tests  will  be  followed  by  simulations  involving  the  New 
England  domain  discussed  in  section  4.  Concurrent  with  this  testing  phase  we  will 
be  developing  the  moisture  parameterizations  allowing  the  production  of  stratiform 
precipitation  in  the  model. 
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In  an  operational  setting,  the  model  would  receive  boundary  conditions  for 
the  CGM  domain  from  one  of  the  National  Weather  Service  operational  models 
(preferably  the  NGM).  These  are  easily  implemented  in  the  current  formulation 
since  the  flow  relaxation  boundary  condition  used  in  the  CGM  allows  direct 
insertion  of  externally  specified  boundary  values.  For  the  initial  real  data  testing, 
we  will  probably  use  objectively  analyzed  observations  taken  during  the  “forecast” 
period  to  create  the  boundary  conditions.  This  should  allow  the  capabilities  of  the 
model  to  be  demonstrated  without  any  bias  due  to  incorrect  boundary  condition 
forecast. 
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APPENDIX 


Derivation  and  Energy  Conservation  of  the  Eta-Coordinate  System 


a.  Derivation  of  the  equations  in  T}~ coordinates 

The  most  common  form  of  terrain  following  coordinate  is  the  cr-system  where 


a 


(A.l) 


where  r  =«  (ps— pt),  and  ps  and  pt  are  the  pressures  at  the  surface  and  the  top  of 
the  model  domain,  respectively.  The  equation  set  in  this  coordinate  system  may  be 
written  (Haltiner  and  Williams  1980) 


^7  +  V*  +  ctolV*  -f  fk  X  V  =  F 

(A.2) 

3d 

_  =  —IT  CL 

(A  .3) 

§2  +  V  ■  TV  +  -  0 

(A. 4) 

Cpdf  _  aw  =  Q 

(A. 5) 

where 

dp 

W  =  dt  =  ra  +  <7'37r/3t  +  v  •  Vir) 

and 

V  =  ui  +  vj  , 

V  =  i3/3x  -f  j3/3y  , 

<t>  =  geopotential, 

a  —  specific  volume, 

/  =  Coriolis  parameter, 

F  =*»  friction  term, 

cP  —  specific  heat  of  dry  air  at  constant  pressure, 

and 

Q  —  heating  term. 
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In  order  to  preserve  conservative  properties  in  the  finite  difference 
equations,  we  write  the  total  derivatives  in  flux  form,  which  can  be  obtained  with 
the  aid  of  the  continuity  equation  (A. 4).  Equations  (A. 2)  and  (A.5)  can  then  be 
written 


3tV  ,  3utV  i  3utV  ,  3irV(7 
3f  '  3z  ^  3y  ^  3ct 


—  ■Kct .itVi  —  / k  X  ■* V  4-  irF 

(A.6) 


3 tT  ,  3 utT  i  BvirT  .  3tT ct  _  7ra  , 

31  r  3i  t  3y  T  3d  Cp  ^  cP 


(A.7) 


The  transformation  to  ^-coordinates  is  carried  out  formally  by  using  the 
definition  of  " 


'-hr  "  = 

so  that 

(j  =  rjH  4-crh 

and  for  some  dependent  variable  A 

fcL4)  =  i  ?Ld 
Verier  H  3 77 


CT  <  CT* 

ct  >  ct* 


(A.8) 


(A  .9) 


f^dl  _ 

1  3  A 

(dsh 

H  377  l 

where  s  —  x,  y,  or  t.  Then  the  ^-coordinate  set  equivalent  to  (A.6),  (A.3),  (A.4),  and 
(A.7)  is 


3 *HV  ,  3utHV  ,  3 vrHW  3 *VHr} 

dt  '  3i  3y  ^  377  “ 

(A. 10) 

-  -  irtfaVoT  -  /k  X  -f-  xHF 


30 

377 


—  Ttfa 


(A. 11) 
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^  +  v.  */ fv  +  i?-0  =  o 

3irHT  ,  auTWT  ,  3utHT  ,  3 *THV  irHa. .  ,  rHQ 
~dT  +  3x  +  +  ~3fj  cTW  +  ~cT 

where  now 

W  =  ^  =  TH77  +  3(CTT)/3f  +  v  •  Vox  (A. 14) 

Note  that  o  is  now  a  dependent  variable  which  is  a  function  of  z,  y,  V  and  t,  and  is 
given  by  (A.9).  As  can  be  seen  by  comparing  (A.10)-(A.13)  with  (A.6),  (A.3),  (A.4), 
and  (A.7),  the  transformation  leads  to  the  prognostic  variables  being  weighted  by  *H 
in  the  T]  system  instead  of  the  tr-weighting  present  in  the  a  system.  Note  also  that 
terms  that  had  been  cr  times  derivatives  of  if  become  derivatives  of  the  quantity 
(cr-ir)  in  the  T]  system. 

b.  Energy  conservation  in  the  r]-aystem 

In  order  to  determine  the  correct  finite  differencing  form  that  will  perserve 
energy  conservation  in  the  finite  difference  equations,  the  energy  conservation 
constraints  of  the  continuous  equations  must  first  be  derived.  The  analysis 
presented  here  closely  follows  that  carried  out  by  Haltiner  and  Williams  (1980)  for 
the  cr-coordinate  equations. 

We  begin  with  the  ^-coordinate  momentum  equation  [equivalent  to  (A.2)l 

^  +  V*  +  aVoT  -!-  /k  X  V  =  F  (A. 15) 

This  is  dotted  with  rHW,  and  rewritten  in  the  flux  form  with  the  aid  of  the 
continuity  equation  (A. 12),  to  yield 

^(ixHK2)  4-  V  •  (1, kHVV 2)  4  ^rHW2)  - 

-irHV  (xHV*  -f  aVair)  +  it HV  •  F  (A. 16) 

The  first  term  on  the  right  hand  side  represents  the  kinetic  energy  production  by 
the  pressure  force.  We  expand  this  term  as  follows 


(A. 12) 

(A. 13) 
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-ttHV  •  (V0  +-  aVo-Tr)  =  -V-  (irHV*)  +  d>V-  (rHV)  -  a tHV  •  Vcrx 


Th  j  continuity  equation  (A. 12)  allows  this  to  be  rewritten 

=  -V  ■  (it HV0)  -  -  aTtfV  •  Vatr 

=  -V  •  (irHV*)  -  -  <xtHV  -  Vctt 

Use  of  the  hydrostatic  equation  (A. 11)  leads  to 

=  -V  •  (*HV<p)  -  <p^P-  -  +  TrHm-Hta)  -  atrtfV  •  Vctt 

Ol  df/ 

Adding  and  subtracting  rHa.(dcrx  /dt),  and  rearranging  yields 
-TTtfV  •  (V*  +  aVCTTT)  =  -V  •  (TTtfV*)  - 


-  H*ap§?  4-  V  •  Vair  +  irHri) 


Then  using  the  definition  of  u,  (A. 14),  and  expanding  the  time  derivatives  of  the 
products  leads  to 


-V  •  (it HV<t>)  - 


3  0tHt? 


-  (H0  -  tf7raa)g  -  *tt^+ 


—  'Ktiouij 

Now,  ( A .  1 1 )  can  be  rewritten  as  /dTl  =  —H{ir era  —  <t>),  and  this  along  with  use 

of  (A. 11)  directly  leads  to 

to  3<J^H77  3<6a  3tt  _3H  „3<f>  3cr 

=  ~v  • (irHV0’  -  3t/ - W  3 1  -  *  IF-  Tai7  at  “ 


Noting  that  ir  is  not  a  function  of  Tj  and  that  3(3cr /3t) /3H  =»  3H  /3t  leads  to  the 
result  that 
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-TtfV  •  (70  +  aVcrir)  =  -V-  (*HV0)  _  ^{<t>*HTl  +  0cr|^) 

(A. 17) 

-  -  ’Hau 

We  form  the  total  energy  equation  by  using  (A. 16)  with  the  substitution 
given  by  (A. 17),  and  adding  it  to  the  thermodynamic  equation  (A. 13)  which  can  be 
written  in  the  form 

&ir HcPT)  +  V-  {ir HVcPT)  +  HcPTHf] )  =  tH(aw  +  Q)  (A. 18) 

di  dr/ 

to  yield 

+  rHcPT )  +  7-  +  *HcPT  -f  t//0) 

+  Mi"'2™  +  TCi’t™  + 

-  ir H(Q  +  V  •  F)  (A. 19) 


If  (A. 19)  is  integrated  from  T)  =  —  1  (where  a  *=  0)  to  T]  =  1  (where  a  -*  1),  we 
obtain 


3. 

3t 


kH(Q  +  V  • 


F)dr? 


(A.20) 


where  we  have  used  Hr)  —  0  at  r)  =»  —1  and  1,  3cr/3f  —  0  at  T)  —  —1  and  1,  and 
0(3ir/3 1)  —  3(0«p»)/3t.  This  is  precisely  the  same  energy  relation  derived  by 
Haltiner  and  Williams  (1980,  their  Eq.  (7-42)J  for  the  ^-coordinate  equations,  if  we 
note  that  HdT)  —  da.  It  is  desirable  for  this  energy  constraint  on  the  continuous 
equations  to  also  hold  for  the  finite  difference  equations.  As  will  be  seen,  this  will 
determine  the  appropriate  form  for  the  finite  difference  equation  set,  and  for  the 
method  of  vertical  finite  differencing. 
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c.  The  energy-conserving  finite  difference  equations 
The  finite  difference  form  of  (A. 10)  is 

|(THVt)  +  ^-iukrHWk)  +  ^<vkrHVk)  + 

^[(Ht/)j.+i/2Vi.  +  l/25  —  (H?7)k-i/2Vfc_i/2)J  =  (A. 21) 

—  rH[y<j>k  +  akV(7feT]  —  /k  X  rHWk  +  ''rHFk 

where  H  takes  on  the  appropriate  value  as  given  by  (A. 8)  depending  on  whether  the 
layer  fc  is  above  or  below  ah.  The  indices  (k— 1/2)  and  (k-fl/2)  refer  to  the  layer 
interfaces  bounding  layer  fc  above  and  below,  respectively.  The  quantity  V  is  the 
interpolation  of  V  from  adjacent  layers  to  the  interface.  Haltiner  and  Williams 
(1980)  show  that  the  advective  terms  will  conserve  both  V  and  (V-V)/2  if 

Vk+i/2  =  ^(Vk+1  +  Vk)  (A. 22) 

We  wish  to  insure  energy  conservation  of  the  rhs  as  well.  The  Coriolis  terms  do 
not  contribute,  and  we  will  ignore  the  friction  term  and  concentrate  on  the  pressure 
gradient  term.  As  with  the  continuous  equations,  we  find  the  rate  of  working  by 
the  pressure  gradient  terms  by  taking  V*-  of  the  first  term  on  the  rhs  of  (A. 21)  to 
obtain 


-rHWk-  [Vtf*  +  akV<7k7r]  =  —  V-  (ir HVk<t>k)  +  <t> ■  (*HVk) 

—<xk-xH\k-  Vctj-T 


(A  .23) 


The  finite  difference  form  of  the  continuity  equation  is 

^  +  V-  rHWk  +  TL[(tf77)k+1/2  -  (Hr»k-l/2]  =  0 


(A. 24) 
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so,  (A.23)  can  be  rewritten 

— *HVk-  [V<£k  -f  akVffkir]  *  — V-  (: wHVK<fiK ) 

—  0k  +  ^[(H?7)k  +  l/2  ~  (W^k-l^J)  —  CL^HV k-  VcTfeT 

Adding  and  subtracting  ^[(Hr?}k+1/20k  +  1/2  —  (H?7)k_i/20k_i/2]  yields 

=  -V-  (7rHVk0k)  -  0k^  -  ^[(Hr7)k+1/20k+1/2  -  (H77)k_l/20k-l/2] 

+  ^[jH7?)k  +  1/2(0k  +  i/2  —  0k)  —  (Hi7)k-l/2(0k-l/2  _  0k)J 

-  akirf/Vk  •  Vcrk* 

Now,  adding  and  subtracting  //•ffa.k(3okn'/3t)  leads  to 

— ir//Vk-  [V*k  +  akVcrkir]  -V-  (ir HVk0k) 

—  ^[(Wi7)k  +  l/20fc+l/'2  —  (W17)k_i/20k_i/2j 

—  (0kW  —  H7rako'k)^  —  0*^^  +  ak-g^- 

(A.25) 

-  +  Vk-  Vcrk7rj 

+  ^[fHi7)k+l),2(0k+1/.2  —  0k)  —  (Hi7)k_1/2(0k_i/2  —  0k)J 
Now,  it  is  easy  to  verify  that 

SP’k  _  O’kS/f  (A  26) 

dt  H  dt 


So,  (A.25)  can  be  written 
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-*WVv-  [V0k  +  a^Vcru-ir]  -  -V-  (rHVk<t>k) 

£^[(^7)*  +  l/20k  +  l/2  (^^?)k  —  l/2^k  —  I/2J 

[o*W  -  H-xa^CT^  -  ^[<t>kH  -  H7rakak]|y 

(A.27) 

-  a.k*H^f  +  Vk  •  Vak7rJ 

2j^£(^^)fc-*-i/2(Ofc+i/2  Ofc)  tH77)k_i/2(^Jc-l/2  —  Ok)] 

But,  using  the  finite  difference  form  of  3(0er)/377  =  — Hi'Ka'OL  —  <f>),  and  noting  that 
3tt /3f  and  3H/3t  are  not  functions  of  T],  we  can  write 


-TtfVk-  [V0k  +  akVCTkTr]  =  -V  -  (TrHVkOk) 


A77[(H»7)k  +  l/20k+l/2  —  (^^7)k  —1  /  20  k— 1/2] 


§*' 
3f 


Ar?[^k  +  ,/2^fc+l/':3F  ~  ^k-l/20fc-l/2 

_ L[m  -^fc-M/23 H  _  1  _p’;c-i/23h1 

Ar?[0k-n/2  w  3(  K-\nr  H  9fJ 

-  +  Vk.  Vffkx] 

^jy£'^^*  +  l/2(^!fc-rl/2  —  Ok)  —  (H77)lc_l/2(0Jc_1 


(A. 28) 


/2 


-  Ok)] 


Using  (A. 26)  again  leads  to  our  final  result 


—  T//Vk-  [V0k  —  akVcrk7r]  =  — V-  (ir/fVk0k) 


_  X^n[^77)k  +  >/20k  +  l/2  —  (WT?)*:  —  i /20fc  —  l /2^ 


-  ^k  +  l/2^+1/2|f  —  <*k-l/  20te  — i/2^yj 

1  fi  „d$k+i/2  -  _3<7k-i/2~| 

-  £^*+1,2*— n  tk-wz*  -- Wt  -  J 

-  +  Vk-  Vo,*] 

+  X^[/^^k+ l/2^k  +  l/2  —  0k)  ~  (H?7)k-l/2(3k-l/2  ~ 


(A. 29) 


On  terra  by  terra  comparison  with  (A. 17),  it  is  clear  that  (A.29)  will  match  the 
continuous  equations,  and  conserve  energy  in  the  same  way  on  summation  over  all 
layers,  if  we  specify  that  U)  is  defined  by 


rHcxkuk  =  ak7 kH  ^  +  Vk-  Vak* 


(A. 30) 


+  !TyI(HV)k  +  l/2{4>k.  +  \/2  —  0k)  —  (ff,7)k-l/2(^k-l/2  —  0k)] 


This  definition  of  w  is  precisely  what  is  needed  to  specify  the  finite  difference 
form  of  the  thermodynamic  and  hydrostatic  equations,  which  we  will  derive  next. 

In  order  to  show  the  constraints  placed  on  the  finite  difference  form  of  the 
thermodynamic  equation,  we  need  to  first  show  two  alternative  finite  difference 
forms  of  the  total  derivitive.  For  a  variable  A  at  level  fc,  we  can  write  the  total 
derivative  in  “flux  form’’  as 


d~KHAk 


T  Vk  •  V 7T H Ak 


+  ^[(H77)k  +  l/2^k  +  l/2  —  (Wi])k_1/2^k-l/2J 


(A.31) 


But  use  of  tiie  continuity  equation  (A. 24)  allows  this  to  be  rewritten  in  an 
“advective  form’’  consistent  with  the  flux  form  as 
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0*  =  Tk/Pk  (A. 35) 

and 

+  Pfc+./2)1/Z/1000j/C  (A. 36) 

The  conservation  of  0  can  be  written  in  advective  form,  (A. 32),  as 

'pa 

rH  ^  4-  Vk  •  V0*  + 

-  9* 

Wf7)fc+i/2(0ic+i /2  ■—  0*.)  +  (Hf])k_l/2{ek  —  0»;-i/2)J  =  o 

Substituting  the  definition  of  0*,  (A. 35),  this  may  be  rewritten 

HI +  Vi  -  +  ’"kslrM +  Vi  -  vhT  (A-37) 

WT7)fc  +  1/2(Pjc0k +•!/;  —  Tk)  +  {HTj)k_l/2(Tk  —  P  *0fc-i/2)J  =  0 

where  we  have  made  use  of  the  fact  that  Pk  can  be  considered  a  function  of  o  and 
•<r.  We  now  add  Ty  times  the  continuity  equation,  (A. 24),  so  that  the  first  term  can 
be  rewritten  in  flux  form,  add  the  finite  difference  form  of  ir3T /dT)  to  both  sides. 
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and  multiply  the  whole  equation  by  cp  to  obtain 


hcpiHT*)  +  V  •  (cP*HTkVk)  +  ^[(Hr7)k+1/2fk+1/2  4  (Hrj)k_l/2f k_l/2)] 


-  'HTT3lfer[i)  +  v‘  '  ^ 


(A.3S) 


4-  -^j\^M‘H)k.+i/2(Ti(,+l/2  —  Pk0k+ 1/2)  4-  (H77)k_1/2(Pk0k_i/2  —  Tk_l/2)J 

The  Ihs  of  (A. 38)  is  the  finite  difference  form  of  the  lhs  of  (A. 18),  and  therefore, 
the  rhs  should  be  equal  to  ttHolOJ  (since  we  are  assuming  Q  =  P  here).  Comparison 
with  (A.37)  shows  that  the  first  terms  will  be  equal  if 


CpTk  3 P k 
Pk  3crk7T 


(A. 39) 


Equating  the  other  terms  leads  to  the  following  relations 


Cp(Tk+i  —  Pk0k+1)  =  0k  — 


Cp(Pkflk_i  —  Tk_i)  =  0k_i  —  #k 


When  (A. 35)  is  used  in  these,  they  can  be  rewritten 


(CpTk+1/2  4-  ^ic+1/2)  —  (CpTk  4-  ^k)  "  P «-.^p(0»l+i/2  —  0k) 


(cPTk  -f-  —  (cpTk_  1/2  4  ^k-1/2)  ”  P fcCP(0k  —  dk  +  1/2) 


(A. 40) 


(A.41) 


We  can  replace  fc  by  Ac  4- 1  in  the  second  of  (A.41),  add  it  to  the  first,  and  again  use 
(A. 35)  to  obtain 


0k  +  l  0*  —  — cP(P  k  +  1  —  P  k)0k 


(A. 42) 


This  represents  a  finite  difference  form  of  the  hydrostatic  equation  that  is 
consistent  with  the  other  equations,  and  provides  e  means  of  calculating  the 
geopotentials  at  all  layers  once  the  geopotential  is  known  at  the  lowest  layer  (where 
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4 


k  =  kbm).  Haltiner  and  Williams  (1980)  show  that  it  is  possible  to  derive  an 
integral  constraint  from  the  energy  conservation  which  yields  the  geopotential  of 
the  lowest  layer  in  terms  of  the  geopotential  differences  of  all  the  other  layers.  It 
is  pointed  out,  however,  that  this  accumulates  the  errors  of  the  layer  calculations 
and  leads  to  large  errors  in  the  geopotential  of  the  lowest  layer.  Experimentation 
with  the  model  verified  this  result.  We  choose  a  simpler  method  of  obtaining  the 
lowest  layer  geopotential  which,  though  not  strictly  consistent  with  energy 
conservation,  yields  accurate  values.  Since  the  boundary  layer  parameterization 
provides  an  “surface-’  temperature,  we  use  this  temperature  to  form  a  “half-layer” 
average  potential  temperature 


_  1/Q 

bm 


+  0 


'l 

s/cJ 


Tnen  an  equation  of  the  form  of  (A. 42)  can  be  used  to  find  in  terms  of  the 

geopotential  at  the  surface. 

The  thermodynamic  equation  can  now  be  written  in  the  form  used  in  the 
model  by  noting  that  (A. 38)  reduces  to 


3  rHTk 
3 1 


+  V  •  7 rtfV*T*  +  ^[(HV)k+l/2Pjk  +  l/7  -  ( HT})k  -I/2P  —  1  /2^j 


*Ho.y  f  3  v 

cp  at  i~  k 


V  (7*77  jr^Qk. 


(A. 43) 


(where  we  have  added  the  diabatic  heating  term).  Thus,  the  equations  (A. 21),  (A .42), 
and  (A.43)  form  a  consistent  finite  difference  set  for  momentum,  geopotential,  and 
temperature  which  conserve  total  energy  in  the  same  manner  as  the  continuous 
equations. 


d.  The  complete  finite  difference  equation  set 

We  restate  here  the  finite  difference  equations  derived  above  and  add  those 
r.ot  yet  discussed  to  form  the  complete  set  used  in  the  model.  The  only  advective 
quantity  in  the  model  whose  finite  difference  form  is  not  discussed  above  is  the  one 
that  governs  specific  humidity,  q.  In  the  absence  of  condensation  or  evaporation, 
however,  specific  humidity  represents  a  conserved  quantity  that  satisfies  dq / dt  «= 
0.  We  write  this  total  derivative  in  the  flux  form  given  by  (A. 31). 
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The  complete  set  of  fimte  difference  equations  for  momentum,  T,  and  q  is 

then 

+  ^(ukirf#Vk)  +  HVk)  + 

^y[(W77)k  +  i/2V)c-(.i/2)  —  (H7?)fc_i/2Vk_i/2)J  =  (A. 21) 

—  irH[V^k  -f-  akVcrkir]  —  /k  X  ir/fVk  4-  *HFk 


0k  +  l  —  4k  =*  — Cp(P  «  +  i  —  P  k)0Jk.  +  l/2 


(A. 42) 


“  +  V‘  '  *]*.*  +  S&*  +  ’"f*  IA'43'> 


and 


where 


— 4-  ^  •  ^HWkQk  +  ^jj\(Hil)k+i/2Qk+i/2  —  (HP)k-i/2<7*-i/zJ  —  F* 

(A. 44) 


0k  m  rK/Pk  (A. 35) 

and 

Pk  -  [i(p*-,/2  4-  pk+l/2)1/2/100oJC  (A. 36) 


Note  that  an  eddy  diffusion  term,  Fk,  is  added  to  the  rhs  of  (A. 43')  and  (A.44)  as 
described  in  section  2  to  help  control  noise.  The  finite  difference  form  of  (2.7), 
which  calculates  the  rate  of  change  of  *  is  written 


3tt 

at 


(A. 45) 


where  H  takes  on  the  correct  values  above  and  below  the  boundary  layer  top  as 
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given  by  (A. 8).  Similarly,  the  finite  difference  form  of  (2.8)  allows  calculation  of 
(Hi 7)  at  each  interface  as 

(Hi 7)*.*i/2  =  — (l?k+ i/2  +  l)  ^ 

(A. 46) 

fc  —  i  *-  -* 

which  can  be  applied  after  (A. 45)  has  been  solved  to  provide  3t/3 t  and  after  the 
boundary  layer  parameterization  has  provided  dh/dt  which  can  be  converted  to 
3<rh/3t.  The  above  finite  difference  equations  are  applied  on  the  staggered  grid 
using  the  averaging  and  differencing  schemes  presented  by  Anthes  and  Warner 
(1978). 
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